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THEORETICAL STUDIES OF A MOLECULAR BEAM GENERATOR 
Introduction 

At present no adequate computer code exists for predicting the effects of thermal 
nonequilibrium on the flow quality of a converging-diverging N 2 nozzle. It is the purpose 
of this research to develop such a code and then perform parametric studies to determine 
the effects of intermolecular forces (high gas pressure) and thermal nonequilibrium (the 
splitting of temperature into a vibrational and rotational-translational excitation) upon 
the flow quality. 

The two models to be compared are given below. See Appendix A for nomenclature 
and additional relationships. 

Model 1 (Equilibrium Model) 

Continuity Equation 


Momentum Equations 


Dp 

Dt 


+ pV • V = 0 


JUV 

pjjf = P v + V * (°ij) 


Energy Equation 


^ + V • ( e t V ) = ^ - W+ Pb • V + V • {oij • V) 
Equation of State 

P=[l~l)pe t , 7=7^’ C P = 

^ v 1 — 1 


( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 
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Model 2 (Nonequilibrium Model) 
Continuity Equation 


Momentum Equations 


Dp - 

Wt +pV - v = 0 


DV t7 _ , x 
Pjfi = pV + V ' ( a u) 


(5) 

( 6 ) 


Energy Equations 

+ V . (e r V) - V • qrt + v • (o,y • F) - 
J* dt (7) 

+ V • (e„?) = - V • q v + pC vv X 

The above equations express the compressible Navier Stokes equations which govern 

the basic dynamics of compressible fluid motion. In these equations we will assume that 

— * 

there are no body forces 6 = 0, and there are no external heat sources Q = 0. In the above 
models the following relations hold 


€t — €f e r — pC xjrfT ~t" pV /2, Cy — pC%j}jTy 


C v — C vr t + C vv \t=T, > C'urt — - R j <f> — 

c ”" = £(£) exp (^)( exp (^) _1 ) 

Qrt = —Krt^T q v = —K v VT v 

K - K + K K -■nr K - 19vk „ - C ' 9cT3/2 

■ti — l\ v + firti K-v — “rt — , > V — 


4m ’ 


T + c 2 


x= T[ 


<j)T 


l-e~*/ T 


. { exp K^ _ ?)‘ 1 ]} 


>-«2(10 n ) (s^Wt) V3 «p [(5^2!i) ,,S _ 
(l+erf(f))sinh(JSaS) 


P(atm)T = 


1/3 


£=(i^!)) ,/6 ( l -(^) 


2/3 N 


( 8 ) 

(9) 

( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 
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where rj is the viscosity given by Sutherland formula with 


ci = 2.16(10 8 ) * 1.488, C 2 = 184.0, and g c = -i2for iST 2 . 
Conservative Equations 

The continuity equation and energy equation are examples of scalar conservation laws 
which have the general form 

r)/v 

(16) 


da — 

— + V -6 = 0 

dt 


where a is a scalar and b is a vector. In Cartesian coordinates we have for the continuity 
equation 

a = p and b = pV x e \ + pV y e 2 + pV z ez (17) 

and for the energy equation we have a = e* and 

3 

b — (g« + P)Vj — VjTxi + q x ei + 


i=i 

3 


( e t d" ■^ > )^2 — ^ 

» = 1 
3 

(et + PWs-^V^ + q, 


i=l 


e 2 + 


®3 


(18) 


where 


e t = pe + p{V ? + V ? + V 3 2 )/2. 

In general orthogonal coordinates (x 1 ,x 2 ,X3) the equation (16) is written in the form 

— ((hih 2 /i3a)) + — — ((^2/1361)) + — — ((^1/1362)) + — — ((hi/1263)) = 0. (19) 

Ot C/X\ C/X2 (sX 3 

where hi,/i 2 ,/i3 are the scale factors obtained from the transformation equations to the 
general orthogonal coordinates. 

The momentum equations are examples of a vector conservation law having the general 

form 


da 

dt 


+ V • (T) = 0 


where a is a vector and T is a symmetric tensor 

3 3 


T — ^ ^ Tjk^j^k. 

k= 1 J=1 


( 20 ) 


( 21 ) 
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In general coordinates we have for the momentum equations 

a = pV and T,y = pV{Vj + P6{j — T{j. 

In general orthogonal coordinates (xi,x 2 ,x 3 ) the conservation law (19) can be written 

r\ O 

— ((/ii/i 2 h 3 a)) + — — ((h 2 /i 3 T • ej)) + — — ((hik s T • e 2 )) + — — ((/ii/i 2 T • e 3 )) = 0 (22) 
at ox i 0x2 0X3 

Example 1 Cartesian Coordinates 

In Cartesian coordinates the model 1 can be written in the strong conservative form 


dU_ dE_ dF_ dG 
dt dx dy dz 


(23) 


where 


E = 


F = 


G = 


P 

pV x 
U = | pVy 
pV z 

et 

pV x 

pVl +P~T XX 
pVxVy ~ T xy 

pV x V z T xz 

(ct 4" P)Vz Vx^xx VyT X y Vz T xz 4“ 9i 
PVy 

pV X Vy T X y 
pVf + P — Tyy 

pV xV z Tyz 

L (e< + P)Vy — V x Ty X ~ Vy T yy ~ ^ z T yz + Qy -I 
PV Z 

rhoV x V z f xz 
pVyVz - T yz 
pV? + P~T ZZ 
L (e t + P)V Z - V X T ZX - VyT zy - V z t zz + q z J 


where the shear stresses are given by 


(24) 


(25) 


(26) 


(27) 


Uj = + Vj ti ) + 6ij\V k ,k 


for i,j,k = 1,2,3. 
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Example 2 Cylindrical Coordinates 

The transformation equations to cylindrical coordinates are given by 

x = r cos 6 
y = r sin 0 

z = z 


and have unit basis vectors 

ei = e r = cos + sin 0e 2 
e 2 = e r = — sin 6ei + cos 0e 2 
^3 = e* = e 3 

and scale factors h\ = 1 , h 2 = r, h$ = 1. Consequently, the continuity equation can be 

|l + « + ^) + M ) =0 (28) 

and the energy equation has the form 

+ p ) Vr ~ VrTrr ~ v$Tre ~ v * Tr * + ^i)) + 

^((r[(e t 4- P)V e - V r T 0r - V B T 9 e ~ V z r 6z + ?*])) + (29) 

^((r[(e t + P)V 2 - V rTzr - V d T z e - V z t zz + q z })) = 0 
The momentum equations have the form 

+ |:((rr„)) + §j{(Tr,)) - Tee + £((r r„)) = 0 
+ §;U'T»r)) + T„ + ~(T„) + §- z ((rT„)) = 0 (30) 

^(MM) + ^Wr.)) + §- $ (T,e) + ^((’■r„)) = 0. 
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For symmetry with respect to the 9 variable we set all derivative terms with respect 
to theta equal to zero and consider only momentum in the r and z directions. Under these 
conditions the only stresses to consider are given by 

r »„ aVr I \/' ia ( rV ') , 

r " + Hr - * - + irj 

r V r 3* / f31 ^ 

** ' dz \r dr dz J 

(dV z dV r \ 

r„=r„ = ^— + — j 

and the above momentum equations reduce to 

§,((rP V r)) + §;(('\P V r + P- ’•rr])) - P + T„ + £((r[pV r V, - T r ,})) = 0 

+ ^((r|pF r ^ - r„|)) + ^((rlpK, 2 + P - r„])) = 0. 

Method of Solution 

In cylindrical coordinates, both the models 1 and 2 can be written in the weak con- 


servative form 


dU dE' dF’ rr „ 

— — -|- —r h ~ — + H — 0 . 

dt dr dz 


The nozzle boundary is described by some defining equation r* = f(z), 0 < z < b and 
therefore we can introduce the transformation equations 


z 

X= b 

r 

y =— = 

r* 


f(xb) 


The computational domain then becomes the x, y domain where 0 < x < 1 and o < y < l 
and the weak conservative form given by equation (33) can be written in the form 

dU dE dF „ 

1T + dI + a^ + H (35) 


where 


dE _ dF' 1 dF_dE“l byf dF' 

dx dx b an dy dy f fdy 


and all other derivative terms, such as in the stresses, are converted to x,y derivatives by 
the chain rule. 
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Using operator slitting we can break the vector equation (35) into a sequence of single 

vector equations which can be used to define the solution. Define the operator L x as 

defining the solution to the vector problem 

dU dE n 

— — + — — — 0 . 
dt dx 

Thus, L x can be defined by the two step predictor-corrector formula 
predictor: Uf] =U,-, - |^(£- +1 , - £*,) 

corrector: U‘- = i [U^ + U% - ^(£?J - £",.,)] 

where Uf] = L X U* -. 

Similarly, we can define an operator L y as the solution of the vector equation 

dU dF n 

~di + Ihj ~ 

where L y is defined as 

predictor: U7] =U'j - 

corrector: 

where Uj* = L y xU ? 

Finally, we define the operator L as the solution of the vector equation 

dU 

m +H = ° 

where L is defined by the modified Euler method or second order Runge-Kutta method 

predictor: U** = — AtH*j 

corrector: U'' = - (/,' - ^-(H^ + tf") 

where U** == L U* . 

We can then define the following sequence of numerical calcuations 

= L x L y LLL y L x U%. (36) 

That is we select a At which satisfies the CFL stability condition 


IAt\ < f IV.I , |V y | 

(Ai) c „ < I _ + — + c, 


(Ar) ! (Ay) 2 


and c = y is the local speed of sound, and then perform the following sequence of 
calculations: 
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Step 1: Solve the system U} t j = L x (7"y . 

Step 2: Solve the system — L y f//y. 

Step 3: Solve the system Ufj = LU 2 j. 

Step 4: Solve the system Ufj — L Uf j . 

Step 5: Solve the system Ufj = L y U* y. 

Step 6: Solve the system I/"* 2 = L x U*j. 

Then redefine 17 "y , go to step 1 and repeat the calculations until the values 17" ■ stop 
changing. This then represents the steady state solution. 

After calculating the vector 

^(steady state) _ [rp, r pV r , r pV zt rct) 


we can calcualte for all i,j values the quantities 


Pi, i 


VriJ 




mi hi 

’■/Pi.i 

mu 

r iPU 


and since 

^( 4 )t J = r j e U,j = r }Pi,3 e i,J + r lPij( V rlj + V *lj)/ 2 
we can determine the internal energy e,y and hence the temperature T,y . In the model 2, 
two temperatures are determined since U will have dimension 5. 

Initial and Boundary Conditions 

We consider a one dimensional nozzle flow where the velocity V z of the gas at any 
cross section is averaged over all cross sectional values. Let A = A(z ) denote the cross 
sectional area of the nozzle. We assume that the mass flow rate across any nozzle section 
is a constant and 


Q o = pV z A = Constant, 


so that we obtain using logaritmic derivatives that 


dp dV, 
— + 

P 


dA 

V, + T=°- 


For isentropic gas flow we let H denote the enthalpy and use the Bernoulli theorem that 

V 2 

— + H = Ho — Constant, 

Zd 
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together with the relations 


, VV 'YJW 

pv' 1 = Constant, v = 1/p, et = C V T = , H = e* + pv = , 

nr - 1 7 - 1 


We can then write 


and use 


dp dp dp 

p dp p 


dp f d P\ _ 2 d P 


7 = f)s=^ — = d/f = -V z dV z 

dp op p 


to obtain 


Using 


7 + f + 


dU* dU 2 


V z v z 


P o 


(-S)-t 


p' h ?i h 

where po and po are the density and pressure while the gas is at rest, we write the enthalpy 
as 


H = 


1 P 


i 

Po_ \ p (7-l )h 


7-1 P 7 - 1 \ Po 
Then the mass rate of flow Qo/A = pV z becomes a function of pressure given by 

' v * = ' |2(ff0 - H)|,/2 “ ^ [tT (cT 2 !)^ 1 - (15) 


From this relation we can plot pV z vs p/po- Note that pV z = 0 when p/po = 1 (gas at rest) 
and p/po = 0 when gas expands into a vacuum. The maximum mass density occurs where 
= 0. This occurs at the throat where V z = c. Since A(z ) is known, we can assume 
an initial value for Qo and then determine two pressures from equation (15) which can be 
plotted as a function of z. We adjust Qo until the two roots are equal at the throat. At this 
value of Qo we can the determine the exit (plenum) pressure p p . If the exit pressure is any 
value greater than this critical value then shock waves can result and the flow could march 
backwards into the nozzle. We use the above analysis to determine starting entrance and 
exit pressures. 
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Our initial conditions axe 
At Entrance 

The temperature, pressure (density), and initial velocities are specified according mass 
flow considerations. The radial velocity is selected in order that the radial velocity in the 
computational coordinates is initially zero. This requires that 

y _ Vjnit y/'(0) Y _ Vjnit 

*° “ \/i + (s//'(0)) 2 " VI + (y/'(0)) 2 ' 

The temperature is initially set to a constant plenum temperature at all interior grid 
points. The plenum pressure is assumed constant at all interior grid points and the radial 
and axial velocities are set to zero. 

Boundary Conditions 

At the entrance (left boundary) the temperature, pressure and velocities axe held 
constant. At the exit (right boundary) the pressure is held as the plenum pressure and the 
temperature and velocities are determined by extrapolation. Along the outside boundary of 
the nozzle the velocities are maintained as zero (no slip boundary condition) and the normal 
pressure and temperatures are maintained as zero. Along the centerline the radial velocity 
is zero while the temperature, axial velocity and density are determined by extrapolation. 
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APPENDIX A 


List of Symbols 


b = Body force vector per unit mass [Newtons/ Kg] 

C vr t = Specific heat at constant volume for rotation-translation [Joule/ Kg K] 

C v = Specific heat at constant volume [Joule/ Kg K] 

C vv = Specific heat at constant volume for vibration [Joule/ Kg K\ 

Dij = Rate of deformation tensor [s -1 ] 

D d -* 

—— = — + V -V Material of substantial derivative 
Dt dt 

e t = Total energy [Joule/m 3 s] 
e, e r ,e v Internal energies [ Joule/m 3 s] 
h — Plank’s constant [ J s] 

Scale factors 


k = Boltzmann’s constant [J/K] 

W 

m — — — Molecular mass [Kg] 

™ a 

N a = Avogadro’s number [mo/ -1 ] 

q r t = Heat input due to rotation and translational energy [ Joule/m 3 s] 
q v = Heat input due to vibrational energy [Joule/m 3 s) 
q = q r t + q v = Total heat energy [ Joule/m 3 s] 

P = Pressure [/Veto/on/m 2 ] 


Ro = Universal gas constant [Joule/ K mole ] 
R = Rq/W = Gas constant [ Joule/ Kg K] 
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List of Symbols 


t = time [5] 

T, T v = Temperatures [ K ] 

V = Velocity [m/s] 

W = Molecular weight of N2 [Kg/ mole] 
r] = Viscosity coefficient [Kg/m s] 

A = Second coefficient of viscosity \Kg/ms\ 

K r t = Thermal conductivity (rotation & translation) [PV/mif] 
K v = Thermal conductivity (vibration) [IV/ms] 
p = Density [Kg/m z \ 
t = ReleLxation time [s] 

T{j = Stresses [Newton /m 2 ] 

4 > = Characteristic vibration [K] 
x,y Computational coordinates 
L x ,L y ,L Operators 
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APPENDIX B 


The nozzle geometry 


NOZZLE GEOMETRY 


is illustrated in the Figure 1 where all dimensions are in millimeters. 



Figure 1. Nozzle Geometry 
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The throat radius is aO = 0.5mm, and the nozzle length is 350.0mm. The left end 
opening is 5.0mm. Starting at the point (0, 5) on the left boundary, we move downward 
on a straight line having a slope of —y/Z. This line intersects a circle of radius 2.0mm 
centered at (a2, 2.5mm) where the point o2 is at the throat and is selected so that the 
slope of the line and the slope of the circle are the same at the point of intersection z — al. 
A point a3 > a2 is selected inorder that a cubic spline can convert the circle into a straight 
line with slope of tan 8°. We write the equation of the circle as 

(r-(a 0 + i?)) 2 + (z-a 2) 2 = R 2 . 

The equation of the line is 

r — 5 = — (tan60°)z. 

At z — a2 we require the slope of the circle to be zero, and at the point al we require that 
the point (al,rl) lie on the circle and further the slope of the circle equals the slope of the 
line. Solving the resulting equations we find that al — a2 = — y/Z, 

al = 5 — (aO + R- v/fl 2 - 3) ^ a2 = ( 8 _( a o + i2_i))/v^. 

For a2 < z < a3 we must find the parameter a3 such that a cubic spline converts the circle 
to a straight line. The cubic spline is represented 

r = S 3 (z) = A{z — a2) 3 -f B{z — a2) 2 + C(z — a2) + D a2 < z < a3 

and is subject to the end conditions that 

S 3 (a2) = aO, Sa(a2) = 0, S3 (a2) = l/R S^aS) = tan 8° Sg (a3) = 0. 


We find that 


A = 


-1 


B = — , C = 0, D = aO. 
2 R' 


12i? 2 tan 8° ’ 

The equation of the nozzle is then given by 

r — f(z) = —\/3 (z — a2) + (aO + R — yJH? — 3 — 3) 

r = f(z) = aO + R — \/.R 2 — (2 — a2) 2 

r = f(z) = S 3 (2) 

r = f(z) = a4 + tan 8 °(z — a3) 


where 


0 < z < al, 
al < z < a2, 
a2 < z < a3, 
a3 < z < a4, 


a4 = S 3 (a3) and a3 = a2 + 2J2tan8°. 
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